#论文中出现的几种多重检验过的代码,输入向量ve和假设个数K，输出拒绝的假设个数，FDR与FNR,1-FDR-FNR
#BH过程
B_H=function(ve,K){
nu=K
rnu=0
snu=0
ve1=sort(ve)
ve2=order(ve)
c2=rep(NA,time=K)
for(i in 1:K){
c2[i]=(i*a)/K}
for(i in K:1){
if(ve1[i]>c2[i]){nu=nu-1}
else{break}}            
if(nu==0){rnu=0
}else{
for(j in 1:nu){
if(Expect[ve2[j]]==0){
rnu=rnu+1
}}}       
m1=K-nu        
if(m1==0){snu=0
}else{
for(j in (nu+1):K){
if(Expect[ve2[j]]!=0){
snu=snu+1
}}}       
if(nu==0){
FDR=0
FNR=snu/K
}
else if(nu==K){
FDR=rnu/K
FNR=0
}
else{
FDR=rnu/nu
FNR=snu/(K-nu)}
return(c(nu,FDR,FNR,1-FDR-FNR))
}

#BL过程
B_L=function(ve,K){
nu=K
rnu=0
snu=0
ve1=sort(ve)
ve2=order(ve)
c2=rep(NA,time=K)
for(i in 1:K){
t=min(1,K/(K-i+1)*a)
c2[i]=1-(1-t)^(1/(K-i+1))}
for(i in 1:K){
if(ve1[i]<=c2[i]){nu=nu-1}
else{break}
}           
if(nu==K){rnu=0
}else{
for(j in 1:(K-nu)){
if(Expect[ve2[j]]==0){
rnu=rnu+1
}}}          
m=K-nu        
if(m==K){snu=0
}else{
for(j in (K-nu+1):K){
if(Expect[ve2[j]]!=0){
snu=snu+1
}}}          
if(nu==0){
FNR=0
FDR=rnu/K
}
else if(nu==K){
FNR=snu/K
FDR=0
}
else{
FNR=snu/nu
FDR=rnu/(K-nu)}
return(c(K-nu,FDR,FNR,1-FDR-FNR))
}

#e-BH过程
e_BH=function(ve,K){
nu=K
rnu=0
snu=0
ve1=sort(ve,decreasing=T)
ve2=order(ve,decreasing=T)
c1=rep(NA,time=K)
for(i in 1:K){
c1[i]=K/(i*a)}
for(i in K:1){
if(ve1[i]<c1[i]){nu=nu-1}
else{break}}         
if(nu==0){rnu=0
}else{
for(j in 1:nu){
if(Expect[ve2[j]]==0){
rnu=rnu+1
}}}        
m=K-nu       
if(m==0){snu=0
}else{ 
for(j in (nu+1):K){
if(Expect[ve2[j]]!=0){
snu=snu+1
}}}          
if(nu==0){
FDR=0
FNR=snu/K
}
else if(nu==K){
FDR=rnu/K
FNR=0
}
else{
FDR=rnu/nu
FNR=snu/(K-nu)}
return(c(nu,FDR,FNR,1-FDR-FNR))
}

#e-down过程
NEW_test=function(ve,K){
nu=K
rnu=0
snu=0
ve1=sort(ve,decreasing=T)
ve2=order(ve,decreasing=T)
c1=rep(NA,time=K)
for(i in 1:K){
c1[i]=K/(i*a)}
for(i in K:1){
if(ve1[i]<c1[i]){nu=nu-1}
else{break}}         
if(nu==0){rnu=0
}else{
for(j in 1:nu){
if(Expect[ve2[j]]==0){
rnu=rnu+1
}}}      
m=K-nu        
if(m==0){snu=0
}else{ 
for(j in (nu+1):K){
if(Expect[ve2[j]]!=0){
snu=snu+1
}}}          
if(nu==0){
FDR=0
FNR=snu/K
}
else if(nu==K){
FDR=rnu/K
FNR=0
}
else{
FDR=rnu/nu
FNR=snu/(K-nu)}
return(c(nu,FDR,FNR,1-FDR-FNR))
}
#p-up-down过程，其中p是阶数，介于1与K之间。
step_up_down=function(ve,K,p){
rnu=0
snu=0
ve1=sort(ve)
ve2=order(ve)
c2=rep(NA,time=K)
for(i in 1:K){
c2[i]=(i*a)/K}
if(ve1[p]>c2[p]){
nu=p
for(i in p:1){
if(ve1[i]>c2[i]){nu=nu-1}
else{break}
}              
if(nu==0){rnu=0
}else{
for(j in 1:nu){
if(Expect[ve2[j]]==0){
rnu=rnu+1
}}}      
m1=K-nu        
for(j in (nu+1):K){
if(Expect[ve2[j]]!=0){
snu=snu+1
}}}
if(ve1[p]<=c2[p]){  
nu=p-1
for(i in p:K){
if(ve1[i]<=c2[i]){nu=nu+1}
else{break}
}           
for(j in 1:nu){
if(Expect[ve2[j]]==0){
rnu=rnu+1
}}     
m1=K-nu       
if(m1==0){snu=0
}else{
for(j in (nu+1):K){
if(Expect[ve2[j]]!=0){
snu=snu+1
}}}  
}
if(nu==0){
FDR=0
FNR=snu/K
}
else if(nu==K){
FDR=rnu/K
FNR=0
}
else{
FDR=rnu/nu
FNR=snu/(K-nu)}
return(c(nu,FDR,FNR,1-FDR-FNR))
}

#p-e-up-down过程,p是阶数，介于1与K之间
e_step_up_down=function(ve,K,p){
rnu=0
snu=0
ve1=sort(ve,decreasing=T)
ve2=order(ve,decreasing=T)
c2=rep(NA,time=K)
for(i in 1:K){
c2[i]=K/(i*a)}
if(ve1[p]<c2[p]){
nu=p
for(i in p:1){
if(ve1[i]<c2[i]){nu=nu-1}
else{break}
}    
if(nu==0){rnu=0
}else{
for(j in 1:nu){
if(Expect[ve2[j]]==0){
rnu=rnu+1
}}}       
m1=K-nu       
for(j in (nu+1):K){
if(Expect[ve2[j]]!=0){
snu=snu+1
}}}
if(ve1[p]>=c2[p]){  
nu=p-1
for(i in p:K){
if(ve1[i]>=c2[i]){nu=nu+1}
else{break}
}        
for(j in 1:nu){
if(Expect[ve2[j]]==0){
rnu=rnu+1
}}     
m1=K-nu    
if(m1==0){snu=0
}else{
for(j in (nu+1):K){
if(Expect[ve2[j]]!=0){
snu=snu+1
}}}        
}
if(nu==0){
FDR=0
FNR=snu/K
}
else if(nu==K){
FDR=rnu/K
FNR=0
}
else{
FDR=rnu/nu
FNR=snu/(K-nu)}
return(c(nu,FDR,FNR,1-FDR-FNR))
}

#通过设置不同的检验水平a，相关系数p，假设个数K与正确假设占比r，随机生成数据并得到检验结果
#表5.1-表5.6数据代码
library(MASSExtra)
library(mvtnorm)
#设置参数
K=1000
r=0.8
n=1
p=0.5
u1=0
u2=3
a=0.05
m=1000*100

Expect=Expect=rep(NA,time=K)
for(i in 1:K){
if(i<=K*r){Expect[i]=u1}
else{Expect[i]=u2}}
  
if(p>-0.5){
var1=matrix(data=p,nrow=K,ncol=K)
var2=(1-p)*diag(nrow=K)
var=var1+var2}
if(p<=-0.5){
var=diag(nrow=K)
for(i in 1:K){
for(j in 1:K){
if(abs(i-j)==1){var[i,j]=p}
}}}

#计算增长因子
b=rep(1,K)
for(i in 1:K){
x=K/i
t1=u2*qnorm(1-a/x,mean=u1,sd=1)-(u2)^2/2
b[i]=x/a/exp(t1)
}
boost1=min(b)
#得到数据表格
ErrorRate=array(data=0,dim=c(m,4,6))
for(i in 1:m){

Sample=mvrnorm(n,Expect,var)
Test_S=Sample
Test_P=Test_S
for(t in 1:K){
if(u2>u1){Test_P[t]=1-pnorm(Test_P[t])
}else{
Test_P[t]=pnorm(Test_P[t])}}

Test_E=Test_S
for(j in 1:K){
Test_E[j]=exp(u2*Test_E[j]-(u2)^2/2)}

Boost_PRDS_E=Test_E*boost1

ErrorRate[i,,1]=B_H(Test_P,K)
ErrorRate[i,,2]=B_L(Test_P,K)
ErrorRate[i,,3]=e_BH(Boost_PRDS_E,K)
ErrorRate[i,,4]=NEW_test(Boost_PRDS_E,K)
ErrorRate[i,,5]=e_BH(Test_E,K)
ErrorRate[i,,6]=NEW_test(Test_E,K)
}

Dmeans=colSums(ErrorRate)/m

Data=data.frame(
BH=Dmeans[,1],
BL=Dmeans[,2],
e-BH=Dmeans[,3],
e-down=Dmeans[,4],
base-e-BH=Dmeans[,5],
base-e-down=Dmeans[,6]
)
row.names(Data)=c("拒绝个数","FDR","FNR","1-FDR-FNR")
Data


#表5.7-表5.12数据代码
library(MASSExtra)
library(mvtnorm)
#设置参数
K=1000
r=0.8
n=1
p=0.5
u1=0
u2=3
a=0.05
m=1000*100

Expect=Expect=rep(NA,time=K)
for(i in 1:K){
if(i<=K*r){Expect[i]=u1}
else{Expect[i]=u2}}
  
if(p>-0.5){
var1=matrix(data=p,nrow=K,ncol=K)
var2=(1-p)*diag(nrow=K)
var=var1+var2}
if(p<=-0.5){
var=diag(nrow=K)
for(i in 1:K){
for(j in 1:K){
if(abs(i-j)==1){var[i,j]=p}
}}}

#计算增长因子
b=rep(1,K)
for(i in 1:K){
x=K/i
t1=u2*qnorm(1-a/x,mean=u1,sd=1)-(u2)^2/2
b[i]=x/a/exp(t1)
}
boost1=min(b)

#得到数据表格Data

ErrorRate=array(data=0,dim=c(m,4,8))
for(i in 1:m){

Sample=mvrnorm(n,Expect,var)
Test_S=Sample
Test_P=Test_S
for(t in 1:K){
if(u2>u1){Test_P[t]=1-pnorm(Test_P[t])
}else{
Test_P[t]=pnorm(Test_P[t])}}

Test_E=Test_S
for(j in 1:K){
Test_E[j]=exp(u2*Test_E[j]-(u2)^2/2)}

Boost_PRDS_E=Test_E*boost1

ErrorRate[i,,1]=step_up_down(Test_P,K,1)
ErrorRate[i,,2]=step_up_down(Test_P,K,K/2)
ErrorRate[i,,3]=step_up_down(Test_P,K,K/4)
ErrorRate[i,,4]=B_H(Test_P,K)
ErrorRate[i,,5]=e_step_up_down(Boost_PRDS_E,K,1)
ErrorRate[i,,6]=e_step_up_down(Boost_PRDS_E,K,K/2)
ErrorRate[i,,7]=e_step_up_down(Boost_PRDS_E,K,3*K/4)
ErrorRate[i,,8]=e_step_up_down(Boost_PRDS_E,K,K)
}

Dmeans=colSums(ErrorRate)/m

Data=data.frame(
K-up-down_1=Dmeans[,1],
K2-up-down=Dmeans[,2],
3K4-up-down=Dmeans[,3],
K-up-down=Dmeans[,4],
1-e-up-down=Dmeans[,5],
K2-e-up-down=Dmeans[,6],
3K4-e-up-down=Dmeans[,7],
K-e-up-down=Dmeans[,8]
)
row.names(Data)=c("拒绝个数","FDR","FNR","1-FDR-FNR")
Data

